
#include <iostream>
#include <cassert>
#include <cmath>
// #include <vector>
#include <utility>
#include <complex>
#include <iomanip>
#include <fstream>
#include <valarray>
using namespace std;

template <typename T>
struct range
{
    const T low, high, step;
    int num;
    range(T l, T h, T st)
        : low{l}, high{h}, step{st}
    {
        num = static_cast<int>((high - low) / step + 1.0);
    }
    valarray<T> operator()()
    {
        valarray<T> val(num);
        int tmp{num - 1};
        for (int i = 0; i < tmp; ++i)
        {
            val[i] = i * step;
        }
        val[tmp] = high;
        return std::move(val);
    }
    int size() const { return num; }
};

int main()
{
    // 一个循环就搞定了  也不需要存储
    constexpr double delta{1e-8};
    constexpr double step{0.1};
    constexpr double high{10.0};
    valarray<double> coordX1(range<double>(0.0, high, step)());
    const int size = coordX1.size();
    // valarray<double> coordX2(range<double>(0.0+delta, 5.0+delta, 1e-3)());
    // sin x^2 的导数
    valarray<double> derivative = 2.0 * coordX1 * cos(coordX1 * coordX1);
    // 前向查分得到的导数
    valarray<double> derivative1 = (sin((coordX1 + delta) * (coordX1 + delta)) - sin(coordX1 * coordX1)) / delta;
    complex<double> cdelta(0, delta);
    cout << cdelta.imag() << endl;
    valarray<complex<double>> coordx2(size);
    for( int i = 0; i < size-1; ++i)
    {
        coordx2[i] = i*step;
    }
    coordx2[size-1] = high;
    valarray<complex<double>> tmp = sin((coordx2 + cdelta) * (coordx2 + cdelta));
    valarray<double> derivative2(size);

    for (int i = 0; i < size; ++i)
    {
        // cout << tmp[i] << "\n";
        // cout << sin(coordx2)
        derivative2[i] = tmp[i].imag() / delta;
    }
    ofstream os{"output.dat"};
    os << "TITLE = \" compare\" \n"
       << "VARIABLES = \"x\", \"theory\", \"forward_diff\", \"complex_deri\"\n";
    os << scientific << setprecision(15);
    for (int i = 0; i < size; ++i)
    {
        os << coordX1[i] << "\t"
           << derivative[i] << "\t"
           << derivative1[i] << "\t"
           << derivative2[i] << "\n";
    }
    os.close();
    cout << sin(complex<double>(0, 1)) << endl;
    return 0;
}